arXiv: 1507.08359vl [math.NA] 30Jul2015 


Structure-preserving integrators for the Benjamin-type equations 

Kimiaki KINUGASA* Yuto MIYATAKE t and Takayasu MATSUO* 

July 31, 2015 


Abstract 

The numerical integration of the Benjamin and Benjamin-Ono equations are considered. They are 
non-local partial differential equations involving the Hilbert transform, and due to this, so far quite few 
structure-preserving integrators have been proposed. In this paper, a new reformulation of the equations 
is stated, and new structure-preserving discretizations are proposed based on it. Numerical experiments 
confirm the effectiveness of the proposed integrators. 


1 Introduction 

In this paper, we consider structure-preserving numerical integration of the Benjamin-type equations. The 
Benjamin equation is represented as 


iif ■ yiix 2 ////-. aLu x fiiixxx — d■ 


where u = ii(t.x), x £ 1, ( > 0. a. /j .y. A are real parameters, the subscript t (or x, respectively) denotes the 
differentiation with respect to the time variable t (or x). The operator L is defined as L = Hd x , where H denotes 
the Hilbert transform 


Hu(x) = -P.V. 
n 


r u{x — y) 
J-oo y 


This equation is often considered on the torus T, i.e. under the periodic boundary condition of length /. In this 
case the Hilbert transform is defined by 


Hu(x ) = yP.V.y cot(^jy S ju(x—y)dy, 


( 1 ) 


or equivalently through its Fourier transform 

F(Hu)(k) = —isgn(&) • F(u)(k). 

The Benjamin equation can formally be seen as a generalization of the KdV equation (a = y = 0) and the 
Benjamin-Ono equation (j8 = y = 0) If5lf28ll. 

The Benjamin equation was first introduced in @ as a governing equation which models unidirectional 
propagation of long internal waves of small amplitude at an interface of two incompressible fluids of different 
density. Global well-posedness is proved for data in both L 2 (M) and L 2 (T) |[24l . Since the KdV and Benjamin- 
Ono equations possess solitary wave solutions of the form u(t,x ) = <p(x — ct), whether the Benjamin equation 
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also has solitary waves has been a major subject. This topic was initially raised by Benjamin f6j |7), and 
extensively studied after that (See, for example, Cl El ED)- Nevertheless, in contrast to the KdV and Benjamin- 
Ono equations, explicit formulae are not known. Since then, numerical studies have been done to approximate 
the solitary waves |[Tl [T5l[T6ll23ll . 

The Benjamin equation has three invariants Il6l 


M = 
1 = 
£ = 


udx 

j 

r i 


-irdx, 


7 2 _k 
2 U 6 1 


t L T CX T l~> 9 \ i 

— —u — — u + — uLu — — u x ) ax. 


p 
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They are also constant when the equation is considered on the torus T. Only the third invariant £ determines a 
Hamiltonian structure u t = d x 8£/8u. 

For PDEs with geometric structures, it is widely accepted that structure-preserving numerical methods often 
yield better numerical solutions than general-purpose methods, especially over a long period of time, and this 
topic has attracting much attention in the last two decades (see e.g. ll22ll for the temporal discretization and iflOl 
ns , 2H |25j for the spatial discretization). For the KdV equation, several structure-preserving methods have 
been proposed; for example, invariants-preserving methods ltl2l[T7li20il . and symplectic and multi-symplectic 
methods £ [4, |30l l. On the other hand, for the Benjamin and Benjamin-Ono equations, only few structure¬ 
preserving methods have been known, because the study on geometry of the Benjamin-type equations are less 
developed, and the non-local operator (H or L) needs to be discretized carefully. We are only aware of a single 
exception: an X-preserving method for the Benjamin-Ono equation, which was proposed by Thomee-Vasudeva 
Murthy ||29l . 

Based on the above observation, in this paper we aim at proposing some new structure-preserving inte¬ 
grators for the Benjamin-type equations. We like to note here that some invariant-preserving discretizations 
are rather obvious in the following sense. First, Thomee-Vasudeva Murthy’s approach mentioned above can 
be readily applied to the Benjamin equation. Second, it is straightforward to derive an ^-preserving inte¬ 
grator based on the Hamiltonian structure mentioned above by utilizing the discrete variational derivative 
method llT2ll20ll2Tll (provided an appropriate discretization of the Hilbert transform, such as the one in Thomee- 
Vasudeva Murthy, or in the present paper). Thus in the present paper we like to try a different approach, which 
is in some sense on top of the literature of the so-called multi-symplectic method. To this end, we first have 
to find a multi-symplectic formulation of the Benjamin type equations. If it is found, it should tell us the local 
behavior of the Benjamin type equations; unfortunately, however, it seems the Hilbert transform (which is a 
non-local operator) prohibits that, at least in a standard manner, and we have to introduce some new ideas. 

In this paper, we do this by extending the concept of the multi-symplecticity so that it can fit into the 
Benjamin-type equations. Then we discretize the equations by the Euler box and Preissmann box schemes. 
Here arises another difficulty—the Preissmann box scheme is generally stabler than the Euler box scheme, and 
thus is more preferable; but as its price it has a disadvantage that it is not uniquely solvable unless the number 
of the spatial grid points is odd (this is caused by an averaging operator in front of the unknown variable). This 
is troublesome in the present context, since in the literature, the discretization of the Hilbert transform has been 
considered only with even number of grid points. In this study, we give a discretization of the Hilbert transform 
also for the odd case, and more importantly, give its theoretical justification. 

This paper is organized as follows. In |Section~2l we briefly review the concept of multi-symplecticity and 
some discretization methods. In [Section 3} we propose new integrators for the Benjamin-type equations. We 


extend the concept of the multi-symplecticity in Section 3TT] discuss the discretization of the Hilbert transform 


and the operator L in |Scction 3.2[ and derive integrators in Sections [373] and [3~4] In Section 4 some numerical 
results are provided. Finally, concluding remarks are given in |Section~5| 

The following notation is used in this paper. For the sake of numerical computation, we impose the periodic 
boundary condition. The domain [0,/] is discretized by uniform meshes with the space mesh size Ax = 1/N. 
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Numerical solutions are denoted by u' n ~ u(iAt.nAx) where At is the time mesh size. When we write only 
the subscript or superscript, it means the associated semi-discretization. We use the abbreviation u n+ \p_ = 
(u n +u n . |_i)/2 (a similar abbreviation is also used for the time index). In order to treat the periodic boundary 
condition, we restrict our consideration to an infinite long vector u = {u n } ne z with the property u n = u n+ #, and 
denote the space to which such periodic vectors belong by X d (in this paper, a bold type always belongs to A c |). 
We use the standard difference operators that approximate d x and d ,: 


8 X u„ = 


8+u' = 


Un -\-1 U n 

Ax 

u i+l -u‘ 
At ’ 


8 t u„ 


8 t u‘ 


ll n Un— 1 

Ax 


At 


8 x u„ — 
8 t u' = 


U„ ■ ! U n —\ 
2Xx ' 
u i+x - u ‘~ 1 
2At 


2 Preliminaries 

In this section, we briefly review the concept of nrulti-symplecticity and some discretization methods. For more 
details on this topic, we refer, for example, to the early references ftO. 271. 


2.1 Multi-symplectic PDEs 

A partial differential equation F(u, u,. u x . u tx ....) = 0 is said to be multi-symplectic if it can be written as a 
system of first-order equations 


Mz t + K Zx = V z S(z), (2) 

with z G a vector of state variables, in which the original variable u is included as one of its components. 
The constant matrices M,K e W 1/d are skew-symmetric, and A is a smooth function depending on z. 

A key property of the nrulti-symplecticity is that there is a multi-symplectic conservation law 

d t co+ d x K = 0, (3) 

where CO and K are differential two-forms defined by 

co = ^dz A Mdz, k = ^dz A Kdz. 

A multi-symplectic PDE also has local conservation laws 

d,E{z) + d x F{z) =0, 
d,l{z) + d x G{z) =0, 


(4) 

(5) 


where 

E(z) = S(z) + l -z]Kz, F(z) = ~\zjKz, 

G(z) = S(z) + ^zJMz, I(z) = -^zjMz. 

Integrating these local conservation laws over the spatial domain, under appropriate boundary conditions and 
appropriate assumptions of F(z) and G(z), leads to the global conservation laws 

s £ = s/ £(z)4l=0 - s I= s/' (2)dv=a <6) 

These quantities are called energy and momentum, respectively. 
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2.2 Multi-symplectic discretizations 


A numerical scheme is called multi-symplectic, if it satisfies a discrete version of the multi-symplectic conser¬ 
vation law As typical multi-symplectic schemes, we give two well-known examples: the Euler box scheme 
and the Preissmann box scheme. 

Let us introduce a splitting of two matrices M and K, i.e. M = M + + M and K = K + + K so that 
M = —M and K \ = —K. The so called Euler box scheme reads 

M+8+4, +M-8 t -i+K + 8+zf n + K-8-z i n = VA(P). (7) 

Although the above splitting is not unique, if we choose M + = \M and K + = \K, the scheme ([7]) is simplified 
to 

M8ti + KS x z i n = V z S(z i n ). (8) 

Hereafter, we consider this special case just for simplicity. The Euler box scheme 0 satisfies the discrete 
multi-symplectic conservation law 


8 , 1 to' + 8 + < = 0. 


where 

K = \k~ 1 A Mdz',, < = U , A Kdzi ■ 

The Preissmann box scheme reads 

M8 t X +k +K8 x + ^=V z S(z* 

The Preissmann box scheme 0 satisfies the discrete multi-symplectic conservation law 

8, + a>‘ n+ i_ +S+Kn + > = 0 , 


(9) 


where 


ffl'i 

;;+j 


= 2 d <+i AMd < + i’ 


K, 


«’+ 

2 dZ " 


1 

2 




1 

2 


The central idea of the Preissmann box scheme is to apply the midpoint rule to the both time and space variables, 
but the Preissmann box scheme differs from the standard midpoint rule in that the domain should be divided 
with odd grid points, i.e. N is odd, for the Preissmann box scheme to ensure the uniqueness of numerical 
solutions. 

In general, multi-symplectic integrators do not inherit the local conservation laws 0, 0 and the global 
conservation laws except for the special cases S(z) is quadratic @. However, it should be noted that 
backward error analysis shows that these global invariants are nearly preserved without any drift ll27l : for the 
Euler box and Preissmann box schemes, errors in energy and momentum are bounded with the order 0(Al 2 + 
Ax 2 ) independently of time. We also note that a semi-discrete scheme may inherit one of the global invariants: 
for example, if we apply the Preissmann box scheme to only the spatial variable, the corresponding semi¬ 
discrete scheme Md t z n +\/i + K8^Zn = V z 5(z n+ j«) inherits the local and global energy conservation laws ll26l . 
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3 Proposed numerical schemes 

3.1 A reformulation of the Benjamin-type equations 

Since the multi-symplecticity tells us both local and global properties of the PDEs, ideally we hope to find 
a multi-symplectic formulation for the Benjamin-type equations. This seems, however, rather demanding, at 
least to the present authors, due to the non-local operators (H and L). Note that the presence of non-local 
operators does not always deny the existence of a multi-symplectic structure. Indeed, some non-local PDEs 
such as the Camassa-Holm and Hunter-Saxton equations possess multi-symplectic structures lfl3l fl4ll . The 
key to the success there is that the non-local operators appear only as some inverse of standard differential 
operators, and can be formally eliminated by multiplying the associated operators; for example, the Camassa- 
Holm equation involves (1 — d 2 )~ l . On the other hand, since the non-locality of the Benjamin equation is due 
to the Hilbert transform, whose inverse is also non-local and thus difficult to treat, it is unlikely that there exists 
a multi-symplectic formulation for the Benjamin equation. 

The above difficulty motivates us to extend the multi-symplectic formulation Q so that the extended for¬ 
mulation can tit into the Benjamin equation. In the formulation Q, V- denotes a standard gradient in the finite 
dimensional setting. We change this gradient with the functional derivative in the infinite dimensional setting, 
and then consider the formulation 



( 10 ) 


where the right hand side denotes the functional derivative of the functional 



with respect to z. This idea is motivated by [0, where a similar approach has been already mentioned. 

Remark 1. More rigorously, in the finite dimensional setting, for S : —> M, the gradient V, in W 1 is defined 

by 


S'(z,rj) = (V z S(z),rj) for all r\ &R d , 


where S' denotes the Gateaux derivative and (■, ■) is the inner product in M rf . In the infinite dimensional setting, 
for S : (L?(T)) d —> M, the gradient, i.e. functional derivative, 8 / 8z in (L 2 (T)) d is defined by 




where S' denotes the Gateaux derivative and (•, ■) is the inner product in ( L 2 (T)) d . For example, if 



the simple calculation 


S(u + Su) — S(u) = f-(uL8u + 8uLu)dx-\-0((8u) 2 ) = I (aLu)8udx +0({8u) 2 ) 



shows that 


—— = aLu. 
bit 
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The Benjamin equation can be written as the formulation (fTO]) with z = [u, v] T , 



0 

1 

2 

0 
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0 
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0 
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_ 1 

0 
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0 

0 

1 

0 

M = 
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, K = 






0 

0 

0 

0 


0 

-1 

0 
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0 

0 

0 

0 


p 

0 

0 
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and 


7 2 * 


S(z) = —wu — 1~U Z — — m 3 + —uLu + ^-v 2 . 

2 6 2 2 

Because of the symmetry of the operator L, the functional derivative is calculated to be 

8S r j , . ,t 

-g— = [—w — /M— ju + aLu, 0, —u, pvj , 

and thus the formulation can be written in the componentwise fashion 

1 A 2 

~(j) t — jSvjc = —w — yu — —u + aLu, 

-~ u t+w x = 0 , 

0 jc — W, 
j8 M.r = j8v. 

Below, we discuss local and global properties of the Benjamin equation based on the formulation (jTOj). 

Let us first consider the multi-symplectic conservation law ([3]). The variational equation associated with 
( jTQ| ) is 

Mdz t + Kdz x = d ^. 

Here, the right hand side is calculated to be 




- 7 -Xu + aL 0—10 
0 0 0 0 
-1 


0 


0 0 0 

0 0/3 


d z- 


Since 


we have 


c Of = ^d zt /\Mdz+ ^dzAMdzt 


Adz - ^dz A ( Kdz x - d. 

2 V V oz 


= — - (dz A Kdz) x — -aLdu A du + -du A aLdu 
= — k x + ad u A Ldu, 


SS 


<*>t + Kx = adn ALdu. 
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( 11 ) 








Here du /\Ldu does not vanish, and thus the multi-symplectic conservation law ([3]) does not hold for the Ben¬ 
jamin equation. However, integrating ( fTT| ) over the spatial domain under the periodic boundary condition, we 
obtain the global property 

d f 

— / (O&x = 0. (12) 

d t J t 

Next we consider the local conservation laws. Taking the inner product of ( fTO] ) with Zt, we obtain 

zjKz x = zJ^ 

because of the skew-symmetry of M. Noticing that 

1 


z Kz x = d x -z Kz - d t \ -z! Kz 


and 


we have 


8S 


a 


a 


z = d t S(z) + —u t Lu - —uLu t , 
oz 2 2 


cc cc 

d t E(z) + d x F(z) = - ^ UtLu + —uLu t . 


Therefore, the local conservation law <0 does not hold. Similarly, we have 


CC cc 

dfl(z) + d x G(z) = ~—u x Lu + —uLu x , 


(13) 


(14) 


and thus the local conservation law Q is not satisfied. Although the local conservation laws do not hold, 
integrating ( fT3] ) and (JT4]) over the spacial domain under the periodic boundary condition immediately leads to 
the global conservation laws 


dt 


dt J t 


dt J t 


7 2 A 


a 


—8 = — E(z)dx=— / --u --u + —uLu — —u x ch = 0, 
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1 


dt dt J t 

Here, the symmetry of the operator L 


—1 = — / 7(z)djc=— / — -u~dLx = 0. 


dt Jj 2 


(15) 

(16) 


/ uLvdx = / (Lu)vdbc 
J T J T 


is used. 

Remark 2. Choosing /3 = y = 0 leads to the formulation (|10 \ for the Benjamin—Ono equation: 
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0 0 0 


u 

-^0 0 
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8z. 
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= — / —wu — —u + —uLu dx 
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Similarly, choosing a = y = 0 leads to the formulation (| 10 )for the KdV equation: 


0^00 
0 0 0 
0 0 0 0 
0 0 0 0 



1 

3 -O- 

i_ 
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w 



V 

t 


0 0 0 -/3 

0 0 10 

0-100 
p 0 0 0 
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w 


V 


8' 

Sz. 


X 3 /3 2 

-wu - —U + — V 

6 2 


cbc, 


hut this coincides with the well-known multi-symplectic form 


since in this case 8S/8z = V z S(z). 
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3.2 Discretizations of the operators H and L 

To solve the Benjamin equation numerically, it is mandatory to discretize the operators H and L = Hd x . In 
particular, the operator L should be discretized so that the symmetry is kept in the discrete setting. 

We first review the approach developed by Thomee-Vasudeva Murthy (29l . However, as will be explained 
soon, their approach only makes sense when the domain is divided into even intervals. Therefore, we shall 
develop a new discretization method for odd intervals so that the Preissmann box scheme is applicable. 

3.2.1 Thomee-Vasudeva Murthy’s approach for H (even intervals) 

In ||29ll a discrete version of the Hilbert transform ([[} for u = {u„} n€ z £ Vi is defined by 



Here, the midpoint rule is used for each interval [ X 2 j,X 2 j+ 2 ] (j = 0,... ,N — 1). This definition is rewritten as a 
discrete convolution 



(17) 


Here we note that {c n } ne z G X L \ and c n = —c_„. Furthermore, we rewrite ( fT7| ) by using the discrete Fourier 
transform. For u = {u n } ne z E W. the discrete Fourier transform F^u is defined by 


{F^u) k = £ M „e- 2 ^, 


12=0 


which also belongs to X^. For v = {vi<}kez £ Vj, the inverse transform v is defined by 



Lemma 1 (l29l). The discrete Hilbert transform defined in © is expressed as 





The following lemma indicates that H^ x u is a second order approximation to Hu. 
Lemma 2 ( ||29ll ). Assume that u is periodic and sufficiently smooth. Then it follows that 


\\H\ x u — Hu\\oo < C(Ar) 2 ||u|| C 3, 


where Hu is an abbreviation of [Hu(xo),Hu(xi),...,Hu(xn~i)\, ||m||oo = max,, \\u„\\ anr/||n|| C 3 = max* \u(x)\ + 
max_ Y | u' (x) | + max v | u" (x) \ + max v | u'" (x) \. 



3.2.2 New approach for H (odd intervals) 


In above, we defined a discrete version of the operator H for even N, and discussed its accuracy. For the Preiss- 
mann box scheme, however, we need its odd number counterpart. Recall that the discrete Hilbert transform is 
written through the discrete Fourier transform in Lemma 1 We wish to obtain a similar expression for odd N. 
For this aim, it is convenient to split (JT|) into two terms: 


Hu(x) = yP.V. 


cot ( y v) u{x — y) dy 


s pv /( co, (s- v )- tan (l- v ))“ (jr -- v,dy 

(L J (t 


= txt lim 

21 e->+0 


(18) 


For the first term, we apply the midpoint rule to each interval [(2 j — 2)Ax,2jAx] (j = I..... (/V — I) /2) and 
the rule f ( ’fix) dx ~ (b — a)f(b ) to the remaining interval [(N — 1)A x,l\. For the second term, we apply 
the midpoint rule to each interval [(2j — l)Ax, (2j + l)Ax] (J = 1 ,...,(N — 1)/2) and the rule f(x) dr ~ 
(.b — a)f(a) to the remaining interval [0, Ax]. We now define a discrete Hilbert transform by 


1 {N ~ 1)/2 / n \ 

(H Ax u) n = — Y cot( — (2j — 1 )Ax)m„ 


y=i 


21 




1 2 Ax 


+ 2/ COt (f ) Un - NAx 

— — tan(0)w„Ax 


1 (A7-l)/2 

~ Yl 7=1 
l (N~ l )/2 

21 

(TV - 1 )/ 2 

21 


Y tan ( Y[ 2jAx) u n 


-27 


2Ar 


Y cot | — (2 / — 
7=1 


~ ^) AX ) Un - 2j+l2AX 
1 / JZ \ 

yj Y tan (^y 2jAx)u n -2j2Ax. 


(19) 


This is rewritten as a discrete convolution 


N -1 


{Haxu)„ = Y dn ~jUj 1 where 

7=0 


^ = f^ c °t(g]), if n is odd, 
l“F tan (w)> if n is even, 


( 20 ) 


and thus it can be expressed through the discrete Fourier transform as desired. 


Lemma 3. The discrete Hilbert transform defined in (20) is expressed as 

HjX\U — F^x ( *Aicld ) h\\M ■ 

where S 0 dd is defined by 

S odd = diag(0, 
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Proof. The proof is similar to that of Lemma 1 Since (H^u) n is written as a convolution in ( [20] ), it immediately 
follows that(F Av // Al «)„ = (F& x d) n (F\ x u) n . Hence, we wish to prove that (F^d)„ = — isgn(n), or equivalently 
FT} (—isgn) = d, where sgn E X d and its components are 


sgn(n) 


'0 if n = 0, 

< 1 if 1 <n<^~, 

-1 if ^ <n<lV-L 
2 — — 


It follows that 

N—\ _■ f(N-\)/2 N -1 \ 

^ 1 (-iig5)„ = -^Sgh(n)e 2 ^ = £ f2nmk/N _ £ & 2nink,N\ 

k= 0 \ *=1 jfc=(JV+l)/2 / 

—2nink/N 


N 2sin(^)cos(^) 

|^ cot (i^)’ if n is odd, 
1 — ^tan(f^), if « is even 


_i (N-D/2 


N 

1 


,2 nink/N 


k= 1 

cos (^) -(- 1 )" 



□ 


Lemma 4. Assume that u is periodic and sufficiently smooth. Then it follows that 

\\HaxU — HuWoc < C(A.r) 2 ||u|| C 3, 


where Huis an abbreviation of [Hu(xo),Hu(xi),... ,Hu(xn~i)\ t , ||«||oo=max„ ||u n || am/||n|| C 3 = max A \u(x)\ + 
max* | u' (x) | + max* | u" (x) | + max* | if" (jc) |. 


Proof. The proof is similar to that of |Lemma 2 
transforms ((p~8|) and (|T9|)) as follows: 


First, we re-express the continuous and discrete Hilbert 


l (Ky\ 

Fhi{x) = — g hm o y cot(^-J (u(x-y) - u(x + y))dy. 

tzj ,o _ 1 .. \ 

21 2 Ax cot ^ 2j J y^n— (2m— 1) ^n+(2m— !))• 


We here define three function \j/(y), and 5(jt,;y) by 


W(y) =ycot(jty), 

u(x-y)-u(x+y) 

Hx,y) = ---, 

8(x,y) = \lf(^<Kx,y) = — cot(^) {u{x-y) - u{x+y)). 

Below, by using these functions, we show that \(H^ x u) n —Hu(x n )\ < C(Ar) 2 ||u|| C 3. Hereafter, we promise that 
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C is a generic constant independently of At and u. First, we rewrite | ( H\ x u) n — Hu{x n )\ as follows. 


| (flxxu)n Hll(x „) | 
! ( N - 1)/2 

21 


Yt L 2Axcot ( ^2/ ‘ ) (“”-(2w-l) ~ M »+(2m-l)) - Yl COt (^7) ( w (•*« “ >’) - « ( x n + >’)) d >’ 


r x 2 m 

rx n 

/ (5(x„,x 2 m -i)-5(x„,y))dy + 

/ (8(x„,x N )-8(x n ,y))dy 

^ x 2m —2 

Jx N -i 


By the standard error estimate for the midpoint and rectangular rules, there exist constants xjm-i < < * 2 m 

(m = 1,..., (N — l)/2) and xn- i < T] < x,v. and we proceed with the estimate 

| (HAxU)n Hli(Xfi) 


M /2 (2At ) 3 5 2 5, e . (Ax ) 2 55, 

JJj 24 5 y 2 (^^-)+ 2 dy * n ’ ^ 


< max 

0 <y</ 


< max 

0<y </ 


5 2 5 


5y 2 

5 2 5 


5y 2 


(•W) 


fej) 


( V )/2 ( A*) 3 

Y + max 

hi 3 


55 


(Av) 3 


/(Ax)* - 

—--h max 

6 0<y</ 


^„,v) 


(Ar) 3 


Here we note that i// is independent of w and thus its norms are bounded. Thus we have 
5 2 5 


max 


5y 2 


(•w) 


<C|| </>(*„, OllcF 


max 

0<y<l 


as , 


0<y </ 

By the Taylor expansion, 

Il0(*»>-)||c* <2 ||m|| C 3, ||0(x„,-)|| c i <2Hc9 

from which we obtain \{H\ x u) n —Hu(x n )\ < C(At) 2 ||»|| C 3 . 


< C\\(j)(Xn,-)\\c l - 


□ 


3.2.3 Symmetric discretization of L 

We now define a discrete version of the operator L = Hd x . Since the discrete Hilbert transform defined above is 
formulated in terms of the discrete Fourier transform, one convenient way for discretizing d x is to employ the 
spectral difference, i.e., for u = { u n } <E X<j, 

DZt* = F Ax l (i^-KjFteU, 

where K € H NxN is defined by 

g = fdiag(0,l,...,f -l,0,-(f — 1),.. -, —1), if N is even, 

\ diag ( 0 , 1 ,... ,F=±, — (A^ 3 -),..., — l), if A is odd. 

Note that this is a real matrix (see, for example, ||19|| ). 

We now define an operator L\ x by 

L Ax u = F Ax ] fhf>s' ] jF AK u. 

where S = S cyc „ if N is even, and S = 5 0 dd if A is odd. This is a symmetric matrix as the next lemma shows, 
which is crucial in the subsequent theoretical analyses. 
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Lemma 5. The operator is symmetric in the sense that L^ x 1 = L\ x . 

Proof. Note that H^ x is a real matrix for both even and odd N, which is clear from the original definition. As 
noted above, the spectral difference operator is also real. Thus the matrix L\ x is real, and it suffices to show 
iftsx)* = Lax, where (• )* denotes the Hermitian conjugate. But it is clear since K and S are both real diagonal 
matrices. □ 


3.3 An Euler box scheme 

We apply the Euler box scheme ([ 8 ]) to the formulation (jTO]) (just for simplicity, we consider the simplified 
version only). To do this, we introduce the notation 


55 

8z 


( 4 ) 


-K - - |K ) 2 + f 4 ( l ax«0«' 

o 

~U l n 

H 


The Euler box scheme for the Benjamin equation reads 

M8 t f n +K8 x j n = ^(z‘ n ), 

We consider if the global property (12] ) is inherited in the discrete setting. 
Theorem 1. For the Euler box scheme it follows that 

8,co' n + 8 x k„ = adu‘ n A (L /Xx du‘) n 


( 21 ) 


and 


N -1 

s+£o>;=o, 

n —0 


where 


K = \k 1 AMdzJ,, < = A Kdzf 


Proof. We consider the discrete variational equation 


M8 t dz l n + K8 x dz l n = d (z l n ) ). 


( 22 ) 


It then follows that 

VK+S'K = d 4 A«W4 + d4A*r4d4 = diAd(f («-)) = «*4 aW^).. 

which immediately indicates ( | 22 | ). □ 

Next we consider the global conservation laws fi3] > and (JT 6 ]). As is the case with standard multi-symplectic 
PDEs and multi-symplectic discretization methods, such invariants cannot be preserved in the fully-discrete 
setting. However, as shown below, semi-discrete schemes possess one of such invariants. 
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Proposition 1. For the semi-discrete scheme 


3S 

MdfZn “I - KSxZ,, — (Zn ) ■ 


it follows that 


N-l 


d t Y, E n = 0, where E„ = S(z „) + z n ) T Kz n . 

n =0 1 


For the semi-discrete scheme 


M8 l z! + Kd x j = ^(j), 


8 t + / Fcbc = 0, where /' = —-{d x z! l ) r Mz' 
J t 2 


it follows that 


Proof We only prove ( |24| ) because ( |25j ) can be proved in a similar manner. 
Taking the inner product ([23]) with d t z n , we obtain 


(■ dtZn) T K8+Zn = ( d t Zn ) ' ^(Zn) 


8S 


Sz 


because of the skew-symmetry of M. Noticing that 


{d,z n ) ] K8+Zn = <5+ ( ^{d t Zn-l) T KZn) ~ d, l '-(8 X ZnVKz, 


and 


we have 


8S 


a 


a 


{df Zn) c (Zn) — 8fS(Z/i) T - {dtUn){E\\U) n U n {L\ x {dtU) ) n , 

8z 2 2 


1 \ Q£ 

dfE n + 8 X ( (d/Zn 1 ) KZnj — {dfUn') (L/SxU) n T — U n {L\ x {d t U) ) n . 


from which ( |24| ) follows. 

Remark 3. For the fidly-discrete scheme the mass conservation 


N-l 

8 ir L “« = 0 

n —0 


is satisfied. 


(23) 


(24) 


(25) 


□ 


3.4 A Preissmann box scheme 

The Preissmann box scheme for the Benjamin equation reads 

MS^p + KStd' 12 = §(£[%)■ 


(26) 
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Theorem 2. For the Preissmann box scheme it follows that 

8 t + <+1/2 + 8+K n +l/2 = adu+Vp A (L At dM , ' +1/2 ), !+ i /2 

and 


where 


$, + L <+1/2 = °> 

n —0 


<+1/2 — 2^ Z " +1 / 2 A ^ Z «+ l / 2 ’ Kn ^ ~ 2 ^ Zn ! ^ ■ 


This theorem is proved in a similar manner to Theorem 1 
Proposition 2. For the semi-discrete scheme 

3S 

Md t Z n+l /2 + K8+Zn = -^(z n + 1 / 2 ), 

it follows that 

N—l 

I 

n—0 

For the semi-discrete scheme 


dt X £ n+l/2 = 0, £„+l/2 = 5 (z„ + i/ 2 ) + - (5+Z„) T ^Z„ + i/ 2 . 


M5+z‘+W' +1/2 = 

oz 


it follows that 


8, + I I'dx = 0, where I‘ = —-(d x z r ) 1 Mz l . 
n 2 


These properties are proved in a similar manner to Proposition 1 


Remark 4. For the fully-discrete scheme (261, the mass conservation 

N -1 

«, + E »:,=<> 

n—0 

is satisfied. 


4 Numerical experiments 

In this section, we check the proposed integrators via several numerical experiments. 

4.1 Solitary wave solution for the Benjamin-Ono equation 

First, we try the Benjamin-Ono equation. We compare the following four integrators. 

• The proposed Euler box scheme ( |2T| ), 

• The proposed Preissmann box scheme ([26]), 
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• The X-preserving scheme by Thomee-Vasudeva Murthy |29l , 

• The Heun method with central difference. 

For the readers’ convenience, the Thomee-Vasudeva Murthy scheme reads 

5 +< + /( 4 + 1 / 2 ) -H^S+S-u? 1 ' 2 = 0 , 

where 

f( u h) = + u h+ u n- l)( M n+l ~ u n- l)- 

Note that this, as well as the Preissmann box scheme, form systems of nonlinear equations at each time step, 
and we need some nonlinear solver. We conducted all the numerical experiments on MATLAB 2010b, and 
employed f solve as our nonlinear solver. The employed Heun scheme reads 

S, + u‘ n = ^ {g{u' n ) +g « + g(u‘ n )At) } , where g(u‘ n ) = -8 X {{u' n ) 2 /2 - LaX} . 

The first three integrators are in some sense structure-preserving, and the last one is not, while keeping the same 
second order accuracy. The last one is employed to see whether actually structure-preserving integrators are 
advantageous. 

We take the parameters to a = I,j8 = 0,y = 0,A = l,l = 30, and At = 2.5 x 10 -3 . We take Ax = 1/256 
for the Thomee-Vasudeva Murthy scheme, which allows only even number of grid points, and Ax = 1/255 
for the rest. The initial data is chosen to u(x n , 0) = 2cA 2 /l — \/\ — A 2 cos (cA(x n — 1/2)), where c = 0.25 and 
A = 2n/cl. The parameter c denotes the speed of the wave. The initial data corresponds to the solitary wave 
solution u(x,t) = 2cA 2 /l — V1 — A 2 cos (cA(x — ct — 1/2)). 

We show the results in Fig. |T]-[3] Fig. [T] shows the wave profiles. The three structure-preserving integrators 
well capture the wave propagation, whereas the Heun scheme exhibits instability with the same discretization 
widths. This confirms the superiority of the structure-preserving integrators. Fig. [2] and Fig. [3] shows the 
evolution of the invariants X and £. From these figures we observe the following two facts. First, the two 
invariants blow up in the Heun scheme, which reflects the instability observed above. Second, the other schemes 
more or less nearly preserve both invariants, at least up to 10 8 . For the Euler box and Preissmann box schemes, 
this is somewhat we expected due to the near conservation properties mentioned above. These properties are 
verified by the backward error analysis (cf. H26JI2ZI). The key of the analysis is that the modified equation is 
also of the form ( flQ| ). The Thomee-Vasudeva Murthy scheme better preserves X, which is a natural result since 
it strictly conserves the invariant by construction. The error observed here should be attributed to the nonlinear 
solver employed in the time-stepping. 

The results above confirms that the proposed integrators are in fact good discretizations. 

4.2 Train of solitary waves in the Benjamin equation 

Next, we try the Benjamin equation. Unfortunately, for this equation, it seems exact solitary wave solutions 
have not been discovered in closed form, whereas many efforts have been done to capture it numerically (see, 
for example, |[Tl[T5l 16l . See also 0.) Thus here we bor row a numerical setting from llT6l Section 2.3], where 
the break down of a Gaussian packet to solitary wave solutions is observed. 

We take the parameters to a = — l,j3 = —1,7 = 1,A = 1,Z = 600, Ax = Z/2048, and At = 10~ 2 . The 
initial data is chosen to u(x n , 0) = 2exp (— (x n — 1 /2) 2 /16). In this experiment, mainly due to the heavy time 
complexity of the Preissmann box scheme, we test only the Euler box scheme. 

The results are shown in Fig. [4] and Fig. [5] Fig. [4] shows the wave profile at t = 100. We observe a train 
of solitary waves, which is exactly what the preceding study lfl6l observed. We here like to note that in lfl6l 
a sophisticated nonlinear scheme was employed, while in the present paper the Euler box scheme is explicit. 


15 




The Euler box scheme 


The Preissmann box scheme 





Figure 1: Solitary wave solution for the Benjamin-Ono equation. 



Figure 2: Evolution of £ in the Benjamin-Ono equation. 
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Figure 3: Evolution of X in the Benjamin-Ono equation. 


Fig. [5] shows the evolution of the invariants. The invariant E is very well preserved. The invariant X deviates 
from the original value in the early phase of wave splitting, and then is well preserved after that. 

Through this experiment, we conclude that the proposed Euler box scheme is in fact a good integrator for 
the Benjamin equation. 


rr(100,x) 



Figure 4: Train of solitary waves in the Benjamin equation. 


4.3 Possible wave breaking in the Benjamin equation 

From the experiment above, we regard the proposed Euler box scheme is a reliable integrator, at least to a 
certain extent. Next, we try to give a new insight about the behavior of the solutions of the Benjamin equation, 
utilizing the integrator. Although in 1241 the global well-posedness of the equation has been proved in X 2 (T), 


17 



















0.06 


o 

fcl 

<u 


0.05 
0.04 
0.03 
0.02 
0.01 
0 

0.01 - 1 - 1 - 1 -- 

0 20 40 60 80 100 



- Energy £ 

-Momentum X 


t 


Figure 5: Evolution of £ and X in the Euler box scheme for the Benjamin equation. 


that does not necessarily prohibits wave breaking in other spaces with higher regularity. The example below 
might suggest such a possibility. 

We take the parameters to a = 0.01, /3 = 0.001 ,7 = 0.1, A = 0.2,1 = 10, Ax = //4096, and At = 10~ 6 . The 
initial data is chosen to u(x n , 0 ) = cos (2kx„/1). 

The results are shown in Fig. [ 6 ]-[ 8 j Fig.[ 6 ]and Fig.|7]show the evolution of u and its derivative u x , respectively. 
In Fig. [ 6 ] it seems the solution develops a steep slope around x ~ 3. Actually in such a place the derivative u x 
seems to blow up (Fig. [7]). Fig. [ 8 ] shows the evolution of the invariants X and £\ we see that both are well 
preserved, which supports that the result is correct. Just to confirm this view, we also tried the fourth order 
Runge-Kutta method and the space discretization with central differences. We observed a similar steep slope 
there (the result omitted here). We also like to point out that in the Runge-Kutta scheme the invariants are not 
preserved up to the same level as the proposed Euler box scheme (Fig.[ 8 ]). This means that even with the fourth 
order temporal accuracy this phenomenon is hard to capture, and the structure-preserving Euler box scheme is 
much advantageous. 



Figure 6 : The solution captured by the Euler box scheme: Evolution of u. 
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5 Concluding remarks 


In this paper, we gave a new reformulation of the Benjamin-type equations, and proposed the Euler box and 
Preissmann box schemes based on the formulation. For the latter scheme, we need the discretization of the 
Hilbert transform on the grids with odd number points, and we provided that with theoretical analysis. The 
numerical experiments confirmed the effectiveness of the proposed structure-preserving schemes. We hope 
that these schemes are useful for better understanding on the behavior of the solutions of the Benjamin-type 
equations. 

One important issue is left unsolved in the present study—although we basically followed the line of the 
discussion in the standard multi-symplectic method, the new formulation is not multi-symplectic in the strict 
sense of the word. In fact, as we saw in Section[3] the standard local conservation laws in usual multi-symplectic 
PDEs are lost, and only their - weaker versions, regarding global invariants integrated over space, are allowed to 
hold. We have to have a deeper understanding about the geometric meaning of the new formulation. 
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